Inputs y paquetes
library(tidyverse)
library(sf)
library(leaflet)
Uruguay <- geouy::load_geouy(c = "Dptos")
iNatUY <- read_csv('../data/observations-175157.csv', guess_max = 25000)
source('../R/funciones_jm.R', encoding = 'UTF-8', local = TRUE)
Para ver completo: descargar versión HTML de este documento y abrir con un navegador.
Este tutorial puede ser instructivo.
#' Columnas iNat a camelCase
#'
#' @param string Ejemplo: 'taxon_family_name'
#'
#' @return string pero en camelCase en vez de separar palabras con '_'
#' @export
#'
#' @examples
#' renombra('scientific_name') # scientificName
#' renombra('taxon_family_name') # taxonFamilyName
renombra <- function(string) {
R.utils::toCamelCase(gsub('_', ' ', string))
}
#' Cantidad de especies nuevas
#'
#' @param nuevas Muestra nueva de especies
#' @param viejas Muestra vieja de especies
#'
#' @return integer. NĆŗmero de especies nuevas
#' @export
#'
#' @examples
#' nuevas_spp(c('a', 'c', 'f'), c('a', 'a', 'c', 'b', 'c'))
nuevas_spp <- function(nuevas, viejas) {
un <- unique(nuevas)
uv <- unique(viejas)
out <- length(un) - sum(un %in% uv)
return(out)
}
((paso pesado))
grid_Uruguay <-
sf::st_make_grid(x = Uruguay, cellsize = 25000, square = F) %>%
st_intersection(.,st_union(Uruguay)) %>%
st_as_sf() %>%
mutate(grid_ID = 1:nrow(.))
last_date <- max(iNatUY$observed_on, na.rm = TRUE)
iNatUY_GIS <- iNatUY %>%
mutate(year = lubridate::year(observed_on),
last_year = observed_on + 365 >= last_date) %>%
filter(captive_cultivated == FALSE,
coordinates_obscured == FALSE) %>%
filter(!is.na(taxon_species_name)) %>%
select(observed_on, year, last_year,
taxon_id,
scientifiName = scientific_name,
class = taxon_class_name,
order = taxon_order_name,
family = taxon_family_name,
genus = taxon_genus_name,
species = taxon_species_name,
decimalLatitude = latitude,
decimalLongitude = longitude,
coordinatePrecision = positional_accuracy,
public_positional_accuracy,
iconic_taxon_name) %>%
sf::st_as_sf(coords = c("decimalLongitude", "decimalLatitude")) %>%
st_set_crs(4326) %>%
st_transform(32721)
Columna last_year: creada para establecer algunos indicadores relacionados a la acumulación de observaciones y registros del último año.
CANTIDAD DE ESPECIES ĆNICAS EN ESA GID
ĆNDICE DE RAREZA DE ESPECIES DE UN ID? DIVERSIDAD ALFA/BETA/GAMA??
ĆNDICE DE REPETICIĆN DE ESPECIES??
REFERENCIA DE LOCALIDADES PARA CADA GID? ((ES DIFĆCIL))
Continuando con la creación del mapa-widget:
# JOIN espacial:
grid_join <- st_join(x = grid_Uruguay,
y = iNatUY_GIS %>%
select(taxon_id, species, class, family,
observed_on, year, last_year,
iconic_taxon_name),
left = TRUE, join = st_contains)
# saveRDS(grid_join, 'data/grid_join.rds')
grid_iNatUY_GIS <-
grid_join %>%
group_by(grid_ID) %>%
# Indicadores:
summarise(
spatialIntensity = ifelse(n_distinct(taxon_id, na.rm = TRUE), n(), 0),
temporalIntensity = n_distinct(year, na.rm = TRUE),
spsList = paste(species, collapse = ';'),
speciesRichness = n_distinct(taxon_id, na.rm = TRUE),
lastYearRecorded = ifelse(spatialIntensity, max(year, na.rm = TRUE), NA),
nObservationsLastYear = sum(last_year),
propObservationsLastYear = nObservationsLastYear / n(),
speciesRichnessLastYear = n_distinct(taxon_id[last_year], na.rm = TRUE),
nNewSpeciesLastYear = nuevas_spp(taxon_id[last_year], taxon_id[!last_year]),
propNewSpeciesLastYear = nNewSpeciesLastYear / speciesRichness) %>%
ungroup() %>%
# Reescalamientos:
mutate(spatialIntensity = ifelse(spatialIntensity,
scales::rescale(spatialIntensity, to = 1:0),
NA),
lastYearRecorded = scales::rescale(lastYearRecorded, to = 0:1),
temporalIntensity = ifelse(temporalIntensity,
scales::rescale(temporalIntensity, to = 1:0),
NA),
indice_prioridad =
scales::rescale(temporalIntensity + spatialIntensity, to = 0:1),
ranking = rank(indice_prioridad,
ties.method = 'max',
na.last = TRUE) %>%
scales::rescale(to = 1:0) %>%
ifelse(is.na(indice_prioridad), NA, .))
# ranking = replace_na(ranking, 1))
# save(grid_iNatUY_GIS, file = 'data/grid_iNatUY_GIS.RData')
Manipular los datos un poco mÔs, y guardar como datos, para simplificar código
datos <-
grid_iNatUY_GIS %>%
select(-spsList) %>%
st_transform(crs = 4326)
# El objeto popup es un vector character con el código HTML
popup <- paste0(
"<strong>GridID: </strong>",
datos$grid_ID,
"<br><strong>Ranking: </strong>",
replace_na(round(100 * datos$ranking, 1), 0), "%",
"<br><strong>Ćndice de prioridad: </strong>",
replace_na(round(datos$indice_prioridad, 2), 1),
"<br><strong>Especies registradas: </strong>",
datos$speciesRichness,
"<br><strong>Especies registradas en el último año: </strong>",
datos$speciesRichnessLastYear,
"<br><strong>Especies nuevas, en el último año: </strong>",
datos$nNewSpeciesLastYear,
" (", round(100 * datos$propNewSpeciesLastYear, 1), " %)",
"<br><strong>Observaciones en el último año: </strong>",
datos$nObservationsLastYear,
" (", round(100 * datos$propObservationsLastYear, 1), " %)"
)
# Proyecto de label para efecto 'mouseover'. Descartado por obstruir la vista:
# mover <- paste0("<strong>Especies registradas: </strong>",
# datos$speciesRichness)
# Objeto (función) para determinar el coloreado de los deptos
# pal <- colorBin(palette = "RdYlGn",
# reverse = FALSE,
# na.color = "#4d0012",
# domain = datos$ranking,
# bins = 8)
# # Otras opciones para la paleta:
# # Esta no me funciona:
pal <- colorQuantile("RdYlGn", datos$ranking, n = 7,
na.color = "#4d0012",
reverse = FALSE)
# # (https://github.com/rstudio/leaflet/issues/94)
# # Parecida, pero con valores continuos:
# pal <- colorNumeric("RdYlGn", reverse = TRUE,
# # c("green", "yellow", "red"),
# datos$ranking)
# El widget de leaflet:
# mapita <-
leaflet(datos) %>%
clearBounds() %>%
addTiles(group = "Open Street Map") %>%
addProviderTiles(providers$OpenTopoMap,
options = providerTileOptions(noWrap = TRUE),
group = 'Open Topo Map') %>%
addProviderTiles(providers$Esri.WorldImagery,
options = providerTileOptions(noWrap = TRUE),
group = 'Imagen') %>%
# addProviderTiles(providers$Stamen.TonerLines,
# # options = providerTileOptions(noWrap = TRUE),
# group = 'Transporte') %>%
# addProviderTiles(providers$Stamen.TonerLabels,
# # options = providerTileOptions(noWrap = TRUE),
# group = 'Transporte') %>%
# # Descartado porque con OSM ya estƔn los departamentos:
# addPolygons(
# data = st_transform(Uruguay, crs = 4326),
# fillColor = NA,
# weight = 2,
# color = 'black',
# fillOpacity = 0,
# group = 'Departamentos',
# ) %>%
addPolygons(
weight = .5,
color = 'white',
fillColor = ~pal(ranking),
fillOpacity = .5,
popup = popup,
# label = mover, # Opción "mouse-over"
highlightOptions = highlightOptions(color = "white",
weight = 5,
fillOpacity = .8,
bringToFront = TRUE),
group = "Grilla"
) %>%
# Control de capas:
addLayersControl(
baseGroups = c("Open Street Map", "Open Topo Map", "Imagen"),
overlayGroups = c("Grilla"), # "Transporte"),
options = layersControlOptions(collapsed = FALSE)
) %>%
# Leyenda de Ćndice de prioridad
addLegend(pal = pal, values = ~ranking, opacity = .5,
na.label = 'Sin registros',
position = 'bottomright', title = 'Ranking', # bins = 8,
group = 'Grilla') %>%
hideGroup("Transporte")
# print(mapita)
# # Obtener datos a partir de 'mapita':
# a <- attr(mapita$x, which = 'leafletData')
# identical(datos, a) # TRUE
# saveRDS(mapita, 'data/mapita.rds')
La curva de acumulación es casera estilo JM, no tiene mucho rigor académico. Para que quede la grÔfica javascript interactiva, es necesario el paquete plotly
p <- grid_join %>%
st_drop_geometry() %>%
arrange(observed_on) %>%
filter(!is.na(taxon_id), !is.na(iconic_taxon_name)) %>%
group_by(iconic_taxon_name) %>%
mutate(c_acum = curva_acum_jm(family),
n_obs = row_number()) %>%
ggplot() +
aes(n_obs, c_acum, color = iconic_taxon_name) +
geom_line() +
ylab('# Familias') +
xlab('# Observaciones')
plotly::ggplotly(p) # Fancy! (interactivo / javascript)